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EXTENSION OF A FINITE VOLUME EXPLICIT TIME MARCHING METHOD TO 
LAMINAR AND TURBULENT FLOW 


ABSTRACT 

This report documents progress made in extending the finite 
volume explicit time marching method to laminar and turbulent 
flow during the time period from January to May 1985. The work 
done is under NASA grant NAG 3-593. Previously, extensions had 
been made to the finite volume method to improve the accuracy of 
the calculation of total pressure in compressible inviscid flow. 
These changes are documented in reference 1 . The current work 
extends these ideas and develops new ideas which allow the 
calculation of laminar and turbulent boundary layers in internal 
flows. The method is verified using four test cases with 
free-stream Mach numbers ranging from .075 to 1.20, 


GOVERNING EQUATIONS 


The unsteady form of the continuity equation, the x-momentum 
equation, and the y-momentum equation, in integral form, are used 
to obtain a steady-state solution through 2-dimensional ducts by 
taking the limit of the unsteady solution as it approaches a 
steady value. The unknown variables are pressure, temperature, 
density, the x-component of the velocity, and the y-component of 
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the velocity. In addition to the continuity and momentum 
equations, we need an equation of state and the energy equation. 
The ideal gas equation of state is used. The energy equation 
currently reduces to the assumption of constant total 
temperature. These five equations are summarized in Table 1. 


CONTROL VOLUMES 

A new control volume has been introduced for this method [1] . To 
eliminate the need for the smoothing of flow properties, there 
must be as many control volumes across the duct as there are 
nodes where these variables are calculated. We need as many 
equations as unknowns. The control volumes also need to be 
located so that errors in continuity and momentum can correctly 
influence the changes in density and velocity without smoothing. 
The current control volume accomplishes this and is shown in Pig. 
1. There are no nodes located along the wall. The nodes are 
located in the middle of the upstream and downstream faces of the 
control volumes. When calculating the flux through a streamwise 
face of an element, the values of the flow properties at the 
node on that face are used. When calculating the flux through a 
cross-stream face, first the properties are calculated at the 
corners of the element using linear interpolation, then the flux 
is calculated using the average of the flow properties at the 
ends of that face. 

The control volumes used by Denton [ 2 ] look like those shown 
in Pig. 2. Fluxes of mass and momentum through each face are 
found by using averages of the flow properties stored at the ends 
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of each face,. However, since the number of nodes (unknowns) is 
greater than the number of control volumes (equations), smoothing 
of flow properties is needed in the crossflow direction to remove 
the dependence of the final solution on the initial guess. 

DISTRIBUTION OF PROPERTIES 

The properties are changed in the flow field after each time step 
because the continuity and momentum equations are not satisfied 
for a given control volume. The amount that properties are 
changed at nodes depends upon the extent to which continuity and 
momentum are not satisfied, the volume of the control volume, and 
the time step. Changes in pressure are distributed to the 
upstream nodes based upon the continuity error. The equation for 
calculating the pressure change is 


(X). 

where $ P is the pressure change, R is the ideal g.as constant, T 
is the local static temperature, and is the density change 
calculated from the continuity error. This density change is not 
used to update the density directly but the density is updated at 
each node using the ideal gas equation of state, 

JL Ca-j )/ETt_ ( 2 ) 

where P T _, is the updated pressure. 
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The lagging correction factor C is defined as, 


Cx-, (3 ) 

and is used to maintain stability as well as to assure that as 
the solution approaches a steady-state the ideal gas equation of 
state will be satisfied at each node, The relaxation factor, P, 
is typically 0.05 and the correction factor is updated every 10 
iterations. It was found that the correction factor must not be 
updated after every step or the solution will become unstable. 
After the pressure is calculated, the momentum equations are 
solved for the changes in velocity. These changes in velocity 
are sent to the downstream nodes. These velocity changes are 
calculated from the momentum error using the updated pressures 
and the old densities and velocities. Finally the static 
temperature is updated using the energy equation. 

When calculating the momentum fluxes, the form of 

the governing equations is approximated by using, 

v-pVV -VCv- p\7) ( 4 ) 

This form subtracts off the continuity error contribution to the 
momentum error and helps the stability characteristics of 
calculations with long thin control volumes. For more details, 
see reference 1, 


4 


BOUNDARY CONDITIONS 


Along the upstream boundary, the total temperature, total 
pressure, and v-velocity are specified for inviscid flow. Along 
the downstream boundary the static pressure is specified. 
Pressures along the solid boundaries are determined from linear 
extrapolation. There is no mass flux across element faces which 
coincide with the solid boundaries. 

For viscous flow, the values of the x-component and 
y-component of velocity are set equal to zero at solid walls. The 
inlet velocity profile is specified along with a freestream total 
pressure. 

INITIAL GUESS 

The initial guess for the inviscid part of the flow field is 
determined from a 1-D inviscid solution. A boundary layer is 
then added along the wall using a constant ratio of boundary 
layer thickness to duct height throughout the duct. The velocity 
profile used in the boundary layer is the inlet velocity profile. 
For certain geometries, an estimate of the blockage effect of the 
boundary layer is used to specify an effective geometry for the 
calculation of the initial solution. 


TIME STEPS 


As outlined in [1], different time steps are used for 
calculating the changes in pressure and velocities. Since we are 
only interested in the steady solution not only can we use 
different time steps at different nodes but also for different 
equations. . If viscous terms are included, the time steps change 
because of the addition of a coefficient from the viscous term. 
The time step is reduced for the calculation of viscous flows. 

The resulting time, or iteration step for the momentum equations 


is 




bum; 

— 


(5) 


and for the continuity equation the iteration step is, 
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In the actual calculations, the time steps are divided by 2 to 4 
because the governing equations are non-linear, 

A simple example can be used to explain in a qualitative 
sense why these different time steps are needed in the 
calculations. The continuity error results in a change in 
pressure at the upstream node associated with a control volume. 
This new pressure is then used in the momentum equation and 
produces a momentum change to improve the continuity error. But 
it can be observed that as the Mach number of the flow decreases 
a given percentage change in pressure produces an increasingly 
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larger change in velocity. The time step used to update the 
pressure must be reduced to keep the changes in velocity stable. 
Table 2 hac a comparison of the ratio of velocity change to 
pressure change with the ratio of time steps calculated using 
equations 5 and 6 for various Mach numbers, Tho quantitative 
agreement of these explains the necessity of different time steps 
for the different equations using the current' updating scheme. 
These numbers are calculated for 1-D compressible flow. 

For low Mach number flows, the time step used in the momentum 
equation can be much greater than that specified by the CFL 
condition. In some of the calculations to be discussed later, 
the time step used in the momentum equation was 500 times the CFL 
condition, while the time step used for continuity was much less 
than the time step specified by the CFL condition. 

CALCULATION OF VISCOUS STRESSES 


The viscous shear stresses are calculated using the following 
equation, 

da. 

(7) 

where yUcW* is the effective viscosity, and du/dy is the local 
velocity gradient. For laminar flow the effective viscosity is 
the absolute viscosity. For the turbulent flow, the viscosity is 
calculated using a Prandtl mixing length model. The equations 
used are summarized in Table 3. Between the near wall point and 
the wall a velocity gradient is calculated and an effective 
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viscosity is calculated at the midpoint using a 

[3], The shear stress at the wall is then calculated from this 
velocity gradient and this effective viscosity. By using this 
effective viscosity, accurate results have been obtained using a 
near wall point at a v+ greater than that usually required . A 
typical value for y+ at the near wall point used in the test 
cases is 30. 

For a non-uniform grid, the top and bottom faces of the 
control volume are not midway between the grid points. The 
velocity gradient calculated from the velocities at the nodes and 
from the distance between these nodes is only valid midway 
between the grid points. In the calculations, the mixing length 
is calculated midway between the points and the shear stress is 
then calculated using this mixing length and the calculated 
velocity gradient. This shear stress is then assigned to the 
face of the control volume between the points. Fig. 3 shows 
schematically the ideas just discussed. This procedure is used 
because the mixing length varies more rapidly than the shear 
stress through the boundary layer . 


NON-UNIFORM GRID SPACING 

For the calculation of internal flows, a numerical algorithm 
must be able to calculate the effect of thin boundary layers 
which grow along the walls of the duct. To be able to calculate 
this type of flow field with a reasonable number of grid points, 
a non-uniform grid must be used with larger spacing in the 
freestream and a higher density of grid points near the walls of 


the duct. Fig, 4 shows how the grid points and control volumes 
might look as you approach the walls of the duct using 
non-uniform grid spacing. The key thing to notice is that the 
aspect ratio (length/height) of the control volumes becomes 
larger as you approach the wall. In the finite-volume 
formulation, this type of control volume requires special 
treatment to maintain the stability of the scheme. 

A more non-uniform grid is needed when calculating a flow 
with turbulent boundary layers since the velocity gradients near 
the wall are much larger. In the present calculations , in the 
boundary layer region, the height of each successive control 
volume is reduced by 50% as the grid approaches the wall. 

Typical distributions of grid points are shown in Table 4 for the 
laminar and turbulent flow calculations. In both cases, the 
boundary layer is only calculated along one wall of the duct. 


A MULTI-VOLUME METHOD FOR PRESSURE CHANGES IN THE BOUNDARY LAYER 

Preliminary calculations of boundary layer flows resulted in 
solutions which became unstable after only a small number of 
iteration steps. After a detailed investigation of the nature of 
this instability, the cause could be attributed to effects 
resulting from the large aspect ratio of the control volumes. The 
first contributing factor was the use of different time steps for 
each control volume and each equation. The time steps that were 
to be used for calculating the change in the density were very 
small. The ratio of the momentum and continuity time steps was 
also very large. The changes in pressure were proceeding in a 


9 


direction which led to large cross velocities and an unstable 
calculation procedure. It was felt that to stabilize this 
calculation procedure, the changes in pressure through the 
boundary layer must be coupled in some manner and that the 
changes in pressure be only dependent on the continuity error and 
not on both the density change through the continuity error and 
the temperature change through the momentum error and its 
resulting velocity change. 

The above realizations resulted in two changes. The first 
change altered the way that the continuity error is used to 
update the flow properties. Previously [1], errors in continuity 
were used to update the density at the node points. The pressure 
was then calculated from the equation of state. An alternative 
procedure has been developed which updates the pressure directly 
from the equation of state (Eqn. 2). To maintain stability, a 
lagging correction factor (Eqn, 3) is used in Eqn. 2 for 
determining the effective pressure to be used in calculating the 
density at a point. In tho limit as the solution reaches steady 
state, the density is evaluated using the correct pressure. This 
procedure is used throughout the flow field and is stable for 
both subsonic and supersonic Mach numbers. 

The second change is to group control volumes in the boundary 
layer to form a larger global control volume. The continuity 
error is calculated for this global control volume and changes in 
pressure are assigned equally to each of the upstream nodes. 
Initially the global control volume extends from the wall' to the 
edge of the boundary layer. Then the global control volume is 
made successively smaller towards the wall. This is shown 
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schematically in Pig. 5. A continuity time step is calculated 
for each global control volume based upon the average properties 
for the control volume. 

TRANSVERSE UPWIND DIFFERENCING 

When the control volumes become long and thin near the wall 
of the duct, the fluxes through the top and bottom faces of the 
control volume become more significant in comparison to the 
fluxes through the streamwise faces. Because the nodes of the 
control volumes are located in the middle of these streamwise 
faces, the predominant flow direction must be in the streamwise 
direction for the discretization method described earlier (see 
Control Volumes) to properly reflect the convective nature of the 
flow?* When the fluxes in the transverse direction become 
significant, the solution procedure may become unstable because 
the diagonal terms in the coefficient matrix become smaller as 
the transverse fluxes increase. This is due to the fact that the 
velocities at the corners of the nodes are determined from 
interpolation. To strengthen the diagonal dominance of the 
coefficient matrix, the momentum fluxes through the transverse 
faces are calculated using the velocities upstream in the 
transverse direction rather than the interpolated values. These 
velocities are multiplied by the mass fluxes through the sides of 
the control volumes to get the total momentum flux. The direction 
of the upwinding is determined from the sign of the continuity 
flux for each face . In the calculations discussed below, this 
upwinding was needed only in the diverging portion of the ducts. 
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TEST CASES 


Pour teat cases have been used so far to test the accuracy 
and stability of the method described in the beginning of this 
progress report. The test cases are of increasing difficulty and 
each serves as a useful check on the method's accuracy. The test 
cases are: 


1. laminar boundary layer in an 
essentially zero pressure gradient. 

2. laminar boundary layer in a favorable 
pressure gradient. 

3. incompressible turbulent boundary layer 
in an adverse pressure gradient. 

4. compressible turbulent boundary layer 
in a transonic diffuser. 

These test cases will now be discussed individually in more 
detail. 

TEST CASE #1 

A laminar boundary layer was calculated in a constant height 
duct. The boundary layer thickness at the inlet was 15% of the 
duct height. The freestream Mach number was 0.43. The inlet 
velocity profile was the Blasius profile. The absolute 
viscosity was 0.01 kg/m s. The duct height was 44 mm and the 
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duct length was 112 mm. The geometry and the grid are shown in 
Pig. 6. The Reynolds number based upon x varies from 5070 to 
11840 along the duct. The duct in 17 inlet boundary layer 
thicknesses long. 

Pig. 7 shows that the development of the velocity profile 
compares very well with that predicted by theory. 


TEST CASE # 2 

A laminar boundary layer was calculated in a converging 
nozzle. The boundary layer thickness was 15% of the inlet duct 
height. The inlet velocity profile was the Blasius profile. 

The inlet height of the duct was 44 mm and the exit height was 
31 mm. The length of the duct was 112mm. The grid and the 
geometry are shown in Pig. 8. The absolute viscosity was 0.01 kg 
/m s. The inlet freestream Mach number is 0.43. 

The calculated velocity profiles were compared with the 
Pohlhausen velocity profiles for the given freestream pressure 
gradients. The velocity profiles are compared at the inlet 
(lambda=0 .0) , midway along the duct (lambda=9 .1) , and the exit 
(lambda=2. 0) . These results are shown in Pigs. 9,10, and 11 
respectively. The agreement is good showing the effect of the 
favorable pressure gradient in changing the shape of the velocity 
profiles ;note that we should only expect qualitative agreement 
with the simple Pohlhausen approximation. 
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TEST CASE # 3 


Incompressible turbulent boundary layer flow in a diverging 
duct was calculated for test case 0141 of the Stanford Conference 
(Samuel and Joubert)[4]. The geometry and grid used by Moore [3] 
are shown in Pig. 12. With this geometry, the top wall is 
treated as inviscid in the calculations. The inlet velocity is 
26 m/s. 

Pig. 13 shows a comparison of the calculated skin friction 
coefficient with the experimental results and with the results 
from the Moore cascade flow program. The agreement is excellent. 
Pig. 14 shows a comparison of the calculated and measured 
velocity profiles at two locations in the duct. The agreement is 
good at x=2.87 m, however, the calculated boundary layer at 
x»3.40 m is thinner than the measurements. It is noteworthy that 
the program has been able to calculate an essentially 
incompressible flow (M<0.075). 

TEST CASE # 4 

Turbulent flow in Sajben's transonic diffuser [5,6] has been 
calculated. The diffuser geometry (Model G) is shown in Pig. 15; 
the throat height was 44 mm and the ratio of the exit height to 
throat height was 1.5. Pig. 15 also shows the computational grid 
which had 83 points in the axial direction and 16 points across 
the flow. The development of a turbulent boundary layer was 
modelled on the curved surface, while the flat wall was treated 
as inviscid. 
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For the calculation, the ratio of the exit static pressure to 
the inlet total pressure was 0.8. This results in transonic flow 
in the diverging portion of the duct with a Mach number of 
approximately 1.235 upstream of a nearly normal shock. In the 
experiment, the flow remained fully-attached throughout the 
diffuser at this test condition. Since the code has not yet been 
developed to handle reverse flow an increased turbulent viscosity 
was used to avoid transient recirculation zones. In particular , 
a fixed minimum value for du/dy (60000 rad/s) was used to 
calculate the turbulent viscosity. The calculation then also 
gave a fully-attached turbulent flow. 

Fig, 16 shows a comparison of the measured and calculated 
wall static pressure distributions on the curved wall. In the 
calculations, the shock is located upstream of the measured 
location, at x/h* = 1.2 compared with x/h* = 1.4. The calculated 
shock location is not sharp and further work is planned to 
improve shock capturing. 

A comparison of the velocity profiles at three axial 
locations along the duct is shown in Fig. 17. The agreement is 
very good, especially downstream of the shock. It may be noted 
that the present calculations give much better agreement with the 
measured velocity profiles than the calculations of reference 6, 
shown as the solid lines. 

Fig. 18 compares the measured and calculated total pressure 
loss for this operating point; the calculation underestimates the 
losses by about 20%. This lower value may be because the 
boundary layer on the flat wall was not included in the 
calculations and may also be due to the slightly weaker shock in 
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the calculation. 


SUMMARY 

The extensions to the explicit finite-volume time marching 
method have made possible the calculation of lam'nar and 
turbulent flows in ducts. The extensions can be summarized as 
follows: 


1. non-uniform grid spacing in the transverse 
direction 

2. a multi-volume method for pressure changes 
in the boundary layer 

3. transverse upwind differencing 

4. different updating procedure for the density 
and the pressure 

5. a Prandtl mixing length model for the turbulent 
shear stresses. 

These extensions have allowed the finite-volume method to 
calculate laminar and turbulent flow for a wide range of 
freestream Mach numbers(0 .075 - 1.20) and geometries. For 
transonic flow in Sajben's diffuser, the calculated overall total 
pressure loss was in reasonable agreement with the measured 
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value. The results suggest that calculations of viscous flow in 
turbomach.tnery blade rows will be possible using the finite 
volume time marching method and that these calculations can model 
the losses. * 
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Continuity 

Ay? - 2^ LfV'dK) At J&V 

it. 

X-momentum ^ 

A(f)Vx) = 2 n (pdAx ^y?Vi V &h/&V 

Y-momentum 

a (^i/y) s ijn /pdA^ + i> Vy v * ) At /ay 

Energy 

To - Corvstcvnt 

Equation of State 

P- pKT 

■#• this is the only viscous term presently modeled 
Table 1 Governing Equations 
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Fig. 2 Typical Control Volumes Used By Denton (2) 
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Table 2 




Ratio of Time Steps for 1-D Compressible Flow 


M 


% velocity change momentum time ste p 

% pressure change continuity time step 


0.02 

1670 

800 

0. 10 

79 

72 

0.80 

1.1 

2.1 

1.00 

0,7 

1.7 
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Table 3 Prandtl Mixing Length Turbulence Model 


- a yf # 

- In the boundary layer 

Jl" smaller of 0.41 y 
or 0.08 8 

h 

- in .41 y region J( a 0.41 y (1- exp^ 

- between the near wall point and the wall 
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Table 4 Typical Transverse Distribution of Grid Points, y/h 


Turbulent Boundary Layer 
§ “ .09 h 

h 13 duct height 


Laminar Boundary Layer 
S m . 15 h 


.0005 

.0020 

.0050 

.0110 

.0230 

.0470 

.0950 

.1910 

.3080 

.4150 

.5210 

.6280 

.7340 

.8410 

.9470 


.015 

.045 

.075 

.105 

.135 

.175 

.225 

.275 

.300 

.334 

.438 

.550 

.650 

.750 

.850 

.950 
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Using Successively Smaller Control Volumes 
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Fig, 10 VELOCITY PROFILE IN CONVERGING NOZZLE 

(LAMBDA - 9.1) 
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13 Skin Friction Coefficient For Samuei and Joubert 



PLOT 4 CASE 0.141 FILES 14.16 


THE 1960-81 AFOSR-dTTM- STAN FORD CONFERENCE ON COMPLEX TURBULENT FLOWS: 
COMPARISON OF COMPUTATION AND EXPERIMENT 
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Pig. 14 Velocity Profiles for Samuel and Joubert 
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Fig. 15 Ceoae try and Grid for Sajben's Diffuser 


CKt >3 

OF POOR QUALnV 


1.0 




0.6 


0.4 



N 


/•* 


cpy 


— $ 


NICHOLSON 


o Experiment 
•— — — Computation 


i.O -zo 


1 
zo 

X/H* 


4.0 


6.0 8.0 




rig. 6 


Surface pressure distribution on upper wall. Model G, 
Rp - «•«. (Too - 10/9, Pr k - P fw 2 - 1.0) 


Fig. 16 Wall Static Pressure Distribution for Sajben's Diffuser 
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Recovery Depends on Shock Strength 



Fig. 18 Total Pressure Loss Comparison For Sajben's Diffuser 


